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ABSTRACT 

A wide variety of astrophysical phenomena involve the flow of turbulent magnetized gas with rela- 
tivistic velocity or energy density. Examples include gamma-ray bursts, active galactic nuclei, pulsars, 
magnetars, micro-quasars, merging neutron stars, X-ray binaries, some supernovae, and the early uni- 
verse. In order to elucidate the basic properties of the relativistic magnetohydrodynamical (RMHD) 
turbulence present in these systems, we present results from numerical simulations of fully devel- 
oped driven turbulence in a relativistically warm, weakly magnetized and mildly compressible ideal 
fluid. We have evolved the RMHD equations for many dynamical times on a uniform grid with 1024 3 
zones using a high order Godunov code. We observe the growth of magnetic energy from a seed field 
through saturation at ~ 1% of the total fluid energy. We compute the power spectrum of velocity 
and density-weighted velocity U = p x ^v and conclude that the inertial scaling is consistent with a 
slope of —5/3. We compute the longitudinal and transverse velocity structure functions of order p up 
to 11, and discuss their possible deviation from the expected scaling for non-relativistic media. We 
also compute the scale-dependent distortion of coherent velocity structures with respect to the local 
magnetic field, finding a weaker scale dependence than is expected for incompressible non-relativistic 
flows with a strong mean field. 
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1. INTRODUCTION 



Turbulence is a fundamental open problem in classi- 
cal physics and is of broad importance in science and 
technology. While non-relativistic turbulence has been 
studied in great detail (e.g., Davidson 2004), relatively 
little attention has been devoted to the properties of tur- 
bulence in relativistic gas. Yet gas flows in astrophysics 
are known to be relativistic and turbulent in a wide va- 
riety of systems under active investigation. Extensive 
evidence of relativistic gas flow now exists, most dra- 
matically in astronomical observations of outflow s from 
gamma-ray bursts (GRBs) (s e e recent reviews by iPiranl 
2004 iFox fc Meszarosl 120061: iWooslev fc Blooml |200€ ' 



Zhang fe Meszarosl 120041: iGehrels et al l 120091) and ac- 
tive galactic nuclei (e.g. Begelman et al.lll984l : IAntonuccll 
[l998TlKrolikl IT9991 : iFrank et all I2002D . Nearer by, rela- 
tivistic outflows are observed from a variety of sources 
including micro-quasa rs, soft gamma repeat ers (SGRs), 
and some supernovae (jSoderberg et al.ll2010ft . Theoreti- 
cal models for the central engines powering a broad range 
of astrophysical phenomena invoke gas accretion onto 
black holes and neutron stars. In addition, the gravi- 
tational wave sources targeted by LIGO involve the hy- 
drodynamical merging of neutron star-neutron star or 
neutron star-black hole binaries. The merging of neu- 
tron star binaries generates shear as the stars touch and 
merge. Shear is also generated as neutron stars are 
shred when merging with a black hole binary companion. 
Kelvin-Helmholtz instability in these shear flows gener- 
ate turbulence which can lead to l arge amplifications of 
magnetic field ()Zhang et al.l I200H) . Such field amplifi- 
cation may be crucial for creating conditions capable of 
extracti ng observed lumino sities from the GRB central 
engine (jGehrels et al.l 



Many astrophysical flows, including the aforemen- 
tioned, are at least partially relativistic (bulk and/or 
thermal Lorentz factor > 1) and all are highly suscep- 
tible to turbulence due to the extremely high Reynolds 
numbers characterizing astrophysical gas. Knowledge of 
the properties of relativistic turbulence is also of general 
importance in physics, with direct applications to early 
universe cosmology, heavy ion colliders, and high energy 
density physics and laboratory plasmas. Research elu- 
cidating the basic properties of relativistic turbulence is 
thus broadly motivated. 

In order to advance knowledge of turbulence in the 
largely unexplored relativistic case, we have begun a se- 
ries of numerical simulations of turbulent flows in mag- 
netized gas with relativistic energy density or velocity. 
In this first in a series of papers on this topic we present 
simulations and analysis of driven turbulence in a trans- 
sonic, super- Alfvenic medium for which the ratio of inter- 
nal to rest mass energy density is of order unity. Differ- 
ences from non-relativistic turbulence may be expected 
since pressure and magnetic fluctuations communicated 
by acoustic and MHD waves modify the fluid inertial 
term in the relativistic case. This is distinct from the 
non-relativistic case where fluid inertia is proportional 
to rest mass density but does not depend on internal 
energy, pressure or magnetic field strength. 

We have studied the same model (SATS1) at a range 
of resolutions up to 1024 3 in order to demonstrate 
numerical convergence. We make comparisons with the 
extensive body of literature available for compressible 
non-rel ativistic hydr o dynam ic and MHD turbulence 
(e.g., iPadoan et al.l 119971; iVazquez-Semadeni et all 
2000; iCho fc Lazarianl 1200% iBeresnvak et all 
Kritsuk et al. 120071 1200a ISchmidt et al.l 



2005 



2008 



2 



J. Zrake and A. I. MacFadyen 



Federrath etaLl 1201131 iLemaster k Stond 120091: 
Burkhart et"a!T l2009t Ikowal k Lazarianl l2010l ). We 
also draw from st udies of Alfvenic turb u lence in the non- 
relativistic (e.g.. iPolitano eTall IT995L ICho k Vishniad 
2000al iMaron k Goldreichl l200ll: iBeresnvak et all 120051: 
Beresnvak k Lazarianl 120091 [2010) and relativistic 
force-free ([Choi 120051 ) regimes. By applying the same 
analyses utilized in these studies to our own models, 
our aim is to elucidate the similarities and differences 
existing between the relativistic and non-relativistic 
limit of MHD turbulence. 

This paper is organized as follows: In <|2]we provide the 
details of our computational scheme and problem setup. 
In Sj3]we report on the saturation of magnetic energy and 
provide power of spectra various quantities of interest. 
We also investigate the scaling of various one and two 
dimensional structure functions of velocity field. In £j4] 
we summarize our findings. 

2. METHODS 

2.1. RMHD formulation 

We have employed Mara, a new unsplit, second-order 
Godunov code which has been written to achieve robust 
and accurate evolution of the RMHD equations on three- 
dimensional rectilinear grids. Mara solves the system of 
ideal RMHD equations in conservation law form. The 
covariant formulation for RMHD can be expressed as a 
system of coupled conservation laws for particle number 
= pu^ and the energy-momentum of the fluid de- 
noted by T^ v = ph*u^u u + p*g» v - WW , where p is the 
rest-mass density of the fluid and u is its four-velocity 
and we use units for which the speed of light c = 1 as 
follows, 



W V N V = 
V„T^ = 

-=Vx(vxB) 



(la) 
(lb) 

(lc) 



Here, W = F^u" is the magnetic field four-vector, and 
h* = 1 + e* + p* /p is the total specific enthalpy, where 
p* = p g + b 2 /2 is the total pressure, p g is the gas pressure 
and e* = e t h + b 2 /2p is the total specific energy density 
with eth being the thermal part. In this study Systcm[T]is 
closed by the adiabatic equation of state, p g = peth (r — 
1) with r = 4/3. Since our equation of state is not 
isothermal are are those used in most previous studies 
(exceptions include [Porter et all ([20021 )1 thermalization 
of the kinetic energy into internal energy results in a 
marginal decline of the sonic Mach number throughout 
the run. This is discussed in more detail in i i3.ll 

The fact that ph* and p* and W directly effect all 
aspects of the energy-momentum dynamics of the fluid 
through their contributions to T^ v means that RMHD 
turbulence dynamics will contain mode couplings not 
present in the non-relativistic case. This is because 
acoustic and RMHD waves will alter eddy dynamics 
via their modifications to fluid inertia, as mentioned 
above. In this study, fluctuations in the thermal en- 
ergy and pressure due to compressive waves will dom- 
inate magnetic fluctuations in T^ v since we consider the 
case with relativistic thermal pressure, pj ' p ~ 1, but sub- 
dominant magnetic field. Alfvenic RMHD turbulence for 



which magnetic energy density is significant or dominant 
(b 2 /2 ~ p) is also of interest and will form the basis for 
a future study. 

The ideal RMHD equations can also be expressed in 
flux conservative form amenable to numerical solution as 



dU <9F J 

dt + ^ ^3 



(2) 



where the conserved quantities 



U = 
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pW 

ph*W 2 -p* - (b ) 2 
ph*W 2 v - b°b 
B 



D 



(3) 



represent the lab frame particle number, total energy (ex- 
cluding rest mass), momentum, and magnetic induction 
respectively, and the corresponding fluxes are 



( Dv* \ 

tv 1 -b°Bi/W +p*vi 
Sv^ - hBi/W + p*V 



(4) 



By introducing the volume averaged conserved quanti- 
ties, 



U 



1 



Xi+l/2 fUj + 1/2 f z k + l/2 



U(x,t)dV (5) 



Zj-l/2 J Vj-l/2 Jz k-l/2 



and rewriting Equation [2] using the divergence theorem, 
we obtain volume averaged time derivatives in each zone, 



(6) 



dt vh3 ' k 
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~ Aa; 


(^i+l/2,j,k 
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~ Ay 


^i,j+l/2,k 






~ Al 







i,j-X/2,k 



where F J , are the fluxes of conserved quantities evaluated 
at the interfaces between cell volumes. Equation [6] would 
in principle complete the integration scheme for advanc- 
ing the solution in time. However, in practice an expres- 
sion which is higher order in time is desirable to maintain 
accuracy. In this study, Mara has been configured to use 
a second order unsplit MUSCL-Hancoc k type integration 
schem e similar to the one described in lMignone k Bodol 
( 2006) , and differs mostly in that magnetic fields are vol- 
ume instead of area-averaged. The scheme is parameter- 
ized around an approximate Reimann solver for obtain- 
ing the intercell fluxes F 3 , and a reconstruction algo- 
rithm for interpolating primitive quantities to the zone 
interfaces. For completeness, the full scheme is described 
below. 

1. Starting with U", invert Equation [3] to obtain P". 
The primitive quantities are P = (p,p g ,v,B) T . 

2. Compute the spatial derivatives in each zone 
using the reconstruction algorithm, and obtain the 
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interpolated primitive quantities at the zone inter- 
faces as follows: 



P J / = P n 



dxi 2 
dP n Ax 3 
dxi 2 



3. Obtain the conserved quantities U^™ +1 ^ 2 and 

U^ 71 " 1 " 1 ^ 2 at the half time-step by applying the 
approximate Riemann solver for the transverse 
fluxes F 2 ' 3 , and a Hancock operator for the nor- 
mal fluxes F 1 . There is one predicted conserved 
state U J ' n+1 / 2 for each direction. For example, in 
the ^-direction, 



U 



1,71+1/2 



At r 



U™ 



At r 



2Ay 
At r. 



2A2 



F 3 (P 



2Ax 

,2,n 

L,i,j+ 

y3,n 



i, fe )~F 2 (P 



2.« 



,2,n 



fl,i,j+l,fc' A L,i,j,k/ 



3,n 



r V r fl,ij,t+1' L,i,i,f!' 



4. Invert Equation[3]to obtain pj>+V 2 from LP v ' l+1 / 2 
for each j = 1,2,3. 

5. Obtain the interpolated primitive quantities at the 
zone interfaces for the half time step. The spatial 
gradients are not recomputed, but are reused from 
the beginning of the time step. 

pi,n+l/2_ p „+l/2 , dP"A^ 

R ~ + dxi 2 

p^n+i/2 =pn+ i/2 dP n Ax 3 
L dxi 2 

6. Complete the time integration by applying the 
fluxes obtained from the Riemann solver in each 
direction. 



TTTl + l _ itn 



At r 

~Ax 
At 

Ay 
At 
~A~z L" 



FVP 



l,n+l/2 pl,n+l/2 



R,i,j,k 

F 2 CP 

r \ r R,i,j,k 



^l-FVP 1 ' 

L,i+l,j,k> \ r R,i+l,j,ki*- L,i,j,k 

2,n+l/2 n A2/p2,n+l/2 p 2,n+l/2x 
L,i,j+l,k) r K r R,i,j+l,ki r L,i,j,k I 

P 3/-p3,n+l/2 p 3,n+l/2 \ _ #3/p3,n+l/2 p3,n+l/2*,' 
r ^R.i.j.fe ' r L,i,j,fe+l-' r \ r R,i,j,k+1' r L,i,j,k I 



1,71+1/2 pl,n+l/2>' 



2.2. Riemann Solvers 



The fluid state in neighboring volumes is formally in- 
terpreted as the initial data of a Riemann problem, 
the solution for which is frequently sought by means 
of an exact or approximate Riemann solver and yields 
the intercell fluxes F 5 . In this study, w e have utilized 
the HLLD approximate Riemann solver (iMignone et al.l 
2009), which computes Godunov fluxes at zone interfaces 
by resolving the fast, Alfven, and contact waves. The 
inclusion of these waves in the Godunov flux approxi- 
mation has been shown to be important for capturing 
the correct concentration of magnetic energy, relative to 



the more diffusive HLLE and HLLC approximate solvers 
(|Beckwith fc Stondl20Tl . 

Due to the fact that Mara employs ILES, the (effec- 
tively numerical) viscous and resistive scales both oc- 
cur near the grid spacing and thus the magnetic Prandtl 
number Pm = Rm/Re ~ 1. However, deviations in 
Pm still occur based on the numerically dissipative be- 
havior of the scheme. Based on qualitative observations, 
the HLLC solver generates roughly the same amount of 
numerical diffusion to the density and pressure fields as 
HLLD, but substantially more to the magnetic field. This 
is interpreted as the HLLD solver achieving a larger nu- 
merical Pm. It is also important to note that as the 
numerical dissipation goes down, the chances of encoun- 
tering robustness issues (see £j2.4p goes up. 

2.3. Magnetic field constraint 

In this study, Mara uses volume-averaged magnetic 
fields which are stored at cell centers. The solenoidal 
constraint V ■ B = (evaluated at cell corners) is main- 
tained to machi ne precision using the constrained trans- 
port method of lTotbl (pOOOl ). 

2.4. Robustness 

Due to the complex nature of the RMHD equations, 
robustness concerns are of great importance. In particu- 
lar, RMHD, and GRMHD codes are known to suffer from 
failures when inverting Equation [3] in order to obtain the 
primitive variables from the conserved ones. This inver- 
sion can fail when the numerical root finder (e.g. secant 
or Newton-Rapheson) fails to locate the root due to an 
insufficiently close initial guess. But it can also fail when 
the root corresponds to a state with negative pressure. 
Overcoming these numerical limitations has been a ma- 
jor source of effort in evolving the models used in this 
study. 

The Mara code has been designed with an extensively 
tested and very robust algorithm for the recovery of prim- 
itive variables. When the failure is strictly numerical in 
nature, it will run through a series of reset values and in- 
version relations. The def ault inversion relat ion used in 
this study is adapted from lNoble et al.l (|2006f ) and in the 
case of an adiabatic equation of state, may be solved us- 
ing a Newton-Rapheson iteration in the single unknown, 
Z = phW 2 . 

2 _S 2 Z 2 + (B-S) 2 (B 2 + 2Z) 



B 2 



{B 2 + ZfZ 2 
^ V) -^2Z^ 



Z -Pg 



The starting value of Z given to the root finder is eval- 
uated from the primitive variables at the previous time 
step. If a suitable solution is not obtained, then the root 
finder is restarted using Z = \JD 2 + S 2 , which becomes 
the exact solution in the limiting case of weak magnetic 
field and small pressure. If a solution is still not ob- 
tained, then the proce dure is repeated usi ng an inversion 
relation adapted from lAnton et al.l (j2006T ) . 



S 2 = (Z + B 



2\2 



w 



1 



D = B 



Z 2 



W 2 
B 2 

2W 2 ' 



- (2Z4 

(B*S t ) 2 
Z 2 
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{B l Sj) 
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Fig. 1. — Slices of the domain taken at t = 12.0T; C , normal to the y-axis. Shown is the rest-mass density (left) scaled linearly between 
0.326 and 1.53, and the (log 10 of) magnetic pressure (right), scaled between 3.16 X 10 -4 and 1.26. 



which can solved using a two-dimensional Newton- 
Rapheson algorithm for the unknowns Z and W . Be- 
cause neither of these solvers contains the other one's 
domain of success, there are circumstances where this 
procedure obtains the solution even when the first solver 
or starting value fails. However, there are states whose 
solution is not obtained by any solver, or whose solution 
is unacceptable due to negative pressure. 

Under these conditions, the code assumes it has inte- 
grated the conserved quantities into an unphysical config- 
uration, and requires further safety procedures in order 
not to crash. The safety feature we have found to be 
the most practical is the addition of small amounts of 
diffusion near "unhealthy" zones. This feature is imple- 
mented as a Lax-Friedrichs flux, where the failed state 
Uij.fe n- U- _ ] k according to 



~U'i,j,k — Ui,j,fc - (f 4 +l/2J,fe 



H-l/2, j.k 



(?) 



where for brevity, we have written the fluxes in the in- 
direction only, and 



H+l/2,j,k 



2d 



(Ui+l,j,k — U 



i,hk)Vi+l/2,j,k 



(8) 



d = 3 is the number of dimensions, and O i+ i^j,k = 1 
if zones (i,j,k) or (i + l,j,k) have been flagged as un- 
healthy, and otherwise. The effect of this prescription 
is to replace U^j ^ with a weighted average of itself and 
the average of its neighboring cells, adding the most dif- 
fusion when r — > 1 and none when r — ¥ 0. Throughout 
this study we have used r = 0.2. This formulation for the 
addition of diffusive terms has several important features. 
Firstly, it naturally obeys the global conservation of U, 
but secondly it obeys the solenoidal magnetic field con- 
straint, because the constrained transport method may 
be applied to the Lax-Friedrichs magnetic field fluxes be- 



fore adding them in Equation [7J 

2.5. Initial conditions and driving 

Our simulations take place in the periodic cube of 
length L centered at the origin. We initialize the do- 
main as a uniform and stationary fluid having rest mass 
density po — 1.0 and gas pressure p g = po/3. We ap- 
ply a uniform magnetic field along the x-direction with 
magnitude 10~ 3 . The flow is driven stochastically on 
large scales according to a prescription we have adapted 
from ISchmidt et all <|2009) • Th e dr iving mechanism is 
intended to mimic the effect of larger flow structures 
in which our domain is embedded, and should thus be 
time-correlated with the turnover time of the largest ed- 
dies, which is of order one light-crossing time of the 
domain, Ti c = L/c. We achieve smooth time corre- 
lation by advancing the Fourier modes a(k, t) of the 
drivi ng field according to an Orn stein-Uhlenbeck pro- 
cess (|Uhlenbeck fc Ornsteinl 119301) <ia(k, t) , which con- 
sists of a restoring force together with a complex- valued 
Gaussian-distributed random- walking term, <iW(k, t): 



da(k,t) =-a(M)^ 

lie 



The projection operator, 




RMS 



<r»(k) 



n 



5p(k) ■ dW(k,t) 



«Py(k)=<V«(k) + (l-C)Vi(k) 



(9) 



is applied to the vector deviate dW(k, t) in order to se- 
lect compressive and vortical driving modes separately, 
according to the parameter £. In this study, we use £ = 1 
which corresponds to a purely vortical driving field. For 
a d etailed study how C ef fects the turbulence statistics, 
see lFederrath et al.l (|2010D . 

The acceleration field is applied to the 4-velocity of 
the flow, at every time step, u(x, t) M- u(x, t) + 
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pended such that e.g., SATS1-1024 refers to the Super- 
Alfenic trans-sonic model 1 at resolution of 1024 3 zones. 




Fig. 2. — (a) Shown are time histories indicating convergence 
of the magnetic energy fraction, the curves represent different res- 
olutions: dotted: 256 3 , dashed: 512 3 , dash-dotted: 768 3 , solid: 
1024 3 . (b) Time histories of various quantities in our highest res- 
olution model, 1024 3 . Shown are the volume- averaged sonic Mach 
number (solid), Alfvenic Mach number (dash-dotted, divided by 
200), and the internal energy (dashed, divided by 5). The Alfvenic 
Mach number at t = YlT lc is ~ 4. 

a(x, t)dt/u°, where the spatial realization is obtained by 
taking the real part of the Fourier mode superposition 



a(x,i) = 5R<( a(k,t)exp(ik-x) } . (10) 

^Q<\U\<K F 

The s pectral profile, a 2 (k) a : k 6 e~ 8k/kl (jVestuto et al.l 
l2003t iLemaster &: Stong 12009ft is normalized to unity 
over the driven wavenumbers. The length scale of max- 
imum driving, l\ — l-Kjkx, and the cutoff If = lit/Kp 
are chosen to be L/4, and L/2 respectively. The small 
subset of driven wavenumbers is chosen for reasons of 
efficiency, since the evaluation of Equation [10] scales 
oc (2Kp + l) 3 and is carried out frequently. The driv- 
ing mechanism which results from this prescription is 
time correlated for TJ C and has RMS power given by 
Prms, which throughout this study has been set to 
0.05, delivering a fractional power per light-crossing 
time, (E tot / E to t}Tic of between 8% and 10% during the 
steady-state period of the run. This relatively mild driv- 
ing power results in a trans-sonic turbulent flow which 
reaches a quasi-steady state after roughly 5 light-crossing 
times. We note that no correlation is enforced between 
real and imaginary parts of c?W(k, t), and thus the driv- 
ing field is statistically helicity-free. We refer to the 
model presented here as SATS1 with the resolution ap- 



3. RESULTS 

3.1. Startup transient and quasi-steady evolution 

Our model of RMHD turbulence begins with spatially 
uniform conditions. Therefore, the early evolution of the 
model is characterized by a startup transient. This tran- 
sient can be roughly partitioned into three stages, which 
have been shown in Figure [2] The forth stage is quasi- 
steady evolution of the model, and is characterized by 
very weak magnetic field growth and a thermalization 
rate balanced by the driving power which causes grad- 
ual increase of the internal energy. The results presented 
in <J3] are obtained from six snapshots taken during this 
phase of the evolution. 

About one correlation time of the driving field (TJ C ) 
passes before its RMS value is reached. During this stage, 
gentle driving gradually increases the mean fluid veloc- 
ity, along with the sonic and Alfvenic Mach numbers. 
For these M ach numbers, we report the relativistic gen- 
eralizations (Ged alinl 1 1993( 1 



M. 



ft 



/3^(1-/3 S 2 A )-V2 
ph 

B 2 



B 2 + 4-n-ph 



(11) 
(12) 

(13) 



where /?/ is the fluid velocity. 

This stage lasts for only about one T; c , until the driving 
field is fully "warmed up." The second stage is charac- 
terized by exponential growth of the magnetic energy, 
and linear growth of kinetic energy. This means that the 
mean Alfven velocity increases more rapidly than the 
bulk velocity, resulting in a sharp drop of the Alfvenic 
Mach number. The transition to quasi-steady evolution 
occurs in the third stage and is complete at 4.7T; C = T sat . 
During the quasi-steady phase, thermalization at small 
scales balances the input of kinetic energy at large scales, 
causing the pressure to gradually increase from 0.333 
to 0.720 over the remaining ~ 7X} C duration of this 
quasi-steady evolution. Also during this stage, magnetic 
structures gather coherence over larger scales, as demon- 
strated in Figure [5jc). The fraction of magnetic to total 
(including thermal) energy increases to 1.5% throughout 
this phase, and the equipartition scale (Pr-(£;) ~ Ps(/c)) 
reaches 1/5 the driving scale. Evolution of a 512 3 model 
(SATS1-512) through 24T /c indicates that the magnetic 
energy fraction stops growing before it reaches 2.0% at 
which time equipartition occurs at 1/3 the driving scale. 
Further discussion of this process is provided in 



3.2. Universality and locality of compressible turbulence 

Universal properties of a turbulent system are those 
which do not dependent upon any boundary, initial, driv- 
ing, or dissipative conditions. The goal of many numeri- 
cal turbulence investigations is to establish universal re- 
lations within turbulent fields at different length scales, 
since by extrapolation they can be used to characterize 
natural systems whose Reynolds numbers are far greater 
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kL/lir 

Fig. 3. — Convergence study for the power spectrum of the veloc- 
ity, Pff(fc). Shown are the averaged power spectra, compensated 
by A: 5 / 3 for the same model at three different resolutions, 512 3 
{dashed), 768 3 {dash- dotted), and 1024 3 {solid). 

than those reached by the numerical model. Universal 
behavior is expected to emerge when the outer and inner 
scales of turbulence are separated far enough that effects 
of driving and dissipation are "forgotten" in an inter- 
mediate range. Exactly what is meant by "far enough" 
depends on the degree with which interactions can take 
place between eddies of disparate sizes. When only ed- 
dies of comparable sizes are able to exchange energy, the 
cascade is said to be highly localized. On the other hand, 
if energy is capable of flowing directly between struc- 
tures separated by large wavenumber, the system is said 
to exhibit nonlocality in its cascade. For such systems, 
enormously large numerical resolution may be required 
in order to resolve universal behavior. 

The locality hy pothesis for incompressible flows 
(|Kolmogorovl 119411 ) led to the well- verified prediction 
that the power in the velocity field, Pff(fc) obeys a power 
law oc fc~ 5 / 3 throughout the inertial range. There are 
now strong arguments that locality holds for compress- 
ible turbulence as well () Aluid [201 If ) . In both compress- 
ible and incompressible MHD, the dynamics of the cas- 
cade are substantiall y more complicated, and nonlocality 
or "diffuse locality" (jBeresnvak fc Lazarianl 120091 120101) 
is expected to hinder the emergence of universal scalings 
in presently available numerical experiments. With re- 
gard to the likelihood that our simulations are capable 
of resolving universality, we s uggest two r e asons to be 
optimistic. The first is that iPorter et al.l ([2002) have 
shown that at A'ls ~ 1 the turbulent energy t rans- 
fer is still consistent with the iKolmogorovl (|1941| ) the- 
ory for incompressible turbulence, meaning that locality 
is likely to be obe yed. The second is th a t the nonlo- 
cality observed by iBeresnvak fe Lazarianl (|2009l) is for 
Alfvenic turbulence in the presence of a strong back- 
ground magnetic field. For our conditions of fully de- 
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kL/2w 

Fig. 4. — Convergence study for the power spectrum of the 
density-weighted velocity, Pjj{k) where U = p 1 / 3 !). Shown are the 
averaged power spectra, compensated by fc 5 / 3 for the same model 
at three different resolutions, 512 3 {dashed), 768 3 {dash- dotted), 
and 1024 3 {solid). 

veloped su per-Alfvenic turbule nce with no substantial 
mean field, iVerma et al.l (|2005l ) predict that the energy 
transfer is local as long as the net helicity is zero. Since 
our driving mechanism is helicity-free, this assumption 
is certainly met. Thus while power law scalings are not 
guaranteed, they should also not be ruled out a priori 
on the basis of distant wavenumber interactions. 

3.3. Convergence study for the power spectrum 

Here we assess the degree to which our simulations 
have resolved the dynamics of the inertial interval by 
studying the power spectrum of the same model at 
different grid resolutions. Careful judgement must be 
used when searching for power law scalings in turbu- 
lence power spectra. Generally, one searches the power 
spectrum for a "flat" interval between two systematic 
over-densities of power at both ends of the spectrum. 
At high wave nu mbers, the so-called "bottleneck" effect 
(|Falkovichlll994f) causes a pile- up of power between the 
inertial and dissipation scales. At low wave numbers, a 
sort of inverse bottleneck occurs as energy rushes away 
from the driving scale. 

We have computed power spectra of the velocity, 
Px(k) and the density- weighted velocity (jKritsuk et al.1 
2007) Pu(k), where U = p 1 /^v for the same model at 
the resolutions N res = 256 3 , 512 2 , 768 3 , and 1024 3 . It is 
important to point out that our interpretation of the res- 
olution study is that the size of the grid spacing is kept 
fixed between runs, so that the box size is reduced to 
L x AV es /1024, and the driving field is moved to higher 
wavenumber. In this interpretation, the power spectra 
are compared to one another at the same number of light- 
crossing times of the largest box size. 

Figure [3] shows the power in the velocity field at each 
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resolution. The main observation is that no inertial 
range is uncovered for the velocity field, but that the 
bottleneck (peaking at kL/2it ~ 10 2 ) is becoming less 
pronounced at higher resolution. Figure |4] on the other 
hand, presents a compelling case that inertial behavior 
has been resolv ed for the densi t y wei ghted velocity p x ^ 3 v 
put fo rward bv iKritsuk et al.l (120071) . and further stud - 
ied in iKowal fc Lazaria nl (12007ft: iSchmidt et al.l (|2008t ); 
iFederrath et al.l (|2010D .~ We believe that this demon- 
strates consistency with the prediction of their simple 
cascade model that Pu(k) oc fc~ 5 / 3 , as well as their nu- 
merical findings for highly c ompressible hydrodynami - 
cal turbulence. We note that iLemaster fc Stonel (|2009l ) 
have also measured Pjj(k) in trans- Alfvenic MHD turbu- 
lence at Ms ~ 6 and obtained a much shallower slope 
of —1.29. We interpret this as evi dence that the cas- 
cade model of IKritsuk et al.l (|2007| ) may be applicable 
to super- Alfvenic, but not trans- Alfvenic flows. This 
interpretation is consi stent with observations made by 
iBoldyrev et al.l (|2002al ) that the low-order hydrodynamic 
statistics of super- Alfvenic turbulence should closely re- 
semble those from the purely hydrodynamical case. 

3.4. Power in compressive versus vortical motions 

As the Mach number of turbulence increases, so does 
the degree of fluid compressibility. This means that more 
of the power in the velocity field is contained in com- 
pressive (dilating or shock-like) structures. In order to 
observe the scale dependence of this effect, we have de- 
composed the velocity field into solenoidal (curl-like) and 
dilatational (divergence-like) parts using a Helmholtz de- 
composition, 



r, 1 2 



P c (k) = |kv k 

p s (k) = |kxv k | 2 



(14a) 
(14b) 



The power Pc{k) in compressive modes is shown in 
Figure [5](a) . Pc(k) follows a power law oc fc -184 
over the wavenumbers k/2ir e [7,31], and no bottle- 
neck is observed. The lack of bottleneck in compres- 
sive m odes has also been observed by (e.g. IPorter et al.l 
[19991) . while a very simila r slope, k 1.79, was observed 
in IFederrath et all pOloh . The ratio P c (k)/P K (k) of 
compressive to total power in the velocity field is shown 
in Figure EJb) . Due to the fact that power is injected 
at large scales using strictly solenoidal modes, there is 
an under-density of compressive power at low wavenum- 
ber. As the details of energy injection are "forgotten" 
at smaller scales, the power in compressive modes grad- 
ually increases toward moderate wavenumber, reaching 
a maximum 9% of the total between the inertial and 
dissipative range. At yet higher wavenumbers through 
the dissipative range, shearing motions become more 
significant, but eventually give way to shocklets, caus- 
ing another rise in the compressive power near the grid 
scale. The overall trend across scales is that ~ 5% of 
the total power is contained in dilata tional motion of the 
fluid, similar to what was found by IPorter et all ([2001 
for Ms ~ 1 flows. It is interesting to note for non- 
relativistic highly supers onic (Ms ~ 10 ) but ot herwise 
very similar conditions, IBoldyrev et al.l (|2002aD found 
that P c (k)/P K (k) ~ 10% - 20% through the inertial 
range, observing a very similar profile to that shown in 




Fig. 5. — Different power spectra taken over the range of times 
Ti c = 8,9, 10, 11, 12 (dash- dotted) and the corresponding time av- 
eraged power spectrum (solid), (a): Pc(k), the power in com- 
pressive velocity modes at wavenumber k, compensated by fc 5 / 3 . 
The blue line has a slope —1.84 and is fit over the wave numbers 
k/2n £ [7,31]. Later times are at lower value than earlier times, 
although the shape of the spectrum is not changing in time, (b): 
Pc(k)/Pi((k), the ratio of compressive to total power in the veloc- 
ity field. Snapshots taken at later times dip lower at moderate to 
high wavenumber. (c): Pg(k)/ Px(k), the ratio of magnetic to ki- 
netic power at wavenumber k. The trend is that this ratio increases 
in time at low and moderate wavenumber. 



Figure Mb). IFederrath et~aT1 (|2010D have studied the ef- 
fect of driving solenoidal versus dilatational modes in 
Ms ~ 5 models the ratios 1/3 and 1/2 respectively. In 
their study, these ratios are constant over the inertial 
range. 

As far as convergence to a quasi-stationary state is 
concerned, we consider the curves in Figure EJa) to be 
robust, because only the overall normalization is chang- 
ing from one snapshot to the next. However, there is a 
definite time trend in the shape of Pc(k)/Px(k), indi- 
cating that further time evolution would be required to 
obtain a time-converged measurement of the compressive 
to total power ratio at different wavenumbers. 

3.5. Power spectrum of magnetic energy 

Figure [5](c) shows the ratio Pb(&)/ PkQz) of power in 
the magnetic field to power in the velocity field. The 
ratio gradually increases from ~ 1% at the driving scale, 
through to super-equipartition at the beginning of the 
dissipative range. It then transitions to become constant 
across scales throughout the dissipative range. Although 
the ratio is still changing in time, with more coherence 
of magnetic structures occurring at large scale, we have 
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Fig. 6. — Shown is the one-dimensional, second order structure 

function of the velocity field S^(£) (x 's) and (+'s). They 

are fit by power laws 8i{£) <x with = 1-03 and rj = 1.02, 
where 21A < i < 82A. 
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Fig. 8. — The slope (~ (boxes) and Qj (diamonds) for differ- 
ent values of the order p. The error bars are obtained by shift- 
ing the fit window by a single bin. Shown also are the pre- 
dictions of jShc & Lcvequc (199 4)) (solid) fo r incom pressible non- 
relativistic hydro t urbule nce, IGrauer et al.l (1993ft (dashed) and 
IMiiller fc BiskamJ (2000 ) (dotted) for incompressible MHD tur- 
bulence. All data have been normalized by £3. 




Fig. 7. — One-dimensional, pth order structure function of the 
velocity field Sp(£), with p = 1 (top) to p = 11 (bottom). 



observed time converged behavior in a 512 3 simulation 
which was run over a longer time. In doing so, we con- 
clude that the interpretation of Figure [5]Jc) is robust. 
After saturation, the ratio of magnetic to kinetic power 
stays constant over the dissipative range, and obeys a 
very straight power law throughout the driving and in- 
ertial ranges, increasing from low to moderate wavenum- 
ber. The scale at which the ratio Pg(fc) / 'Pjf(fc) = 1 
occurs at k ~ 5Kp (see H2.5[) in the most recent snap- 
shot. However, the trend is for the equipartition scale to 
move to lower wave number as magnetic structu r es form 
coherence over larger scales. iCho fc Vishniad (|2000bh 
have observed that once fully steady state is reached, 
the equipartition scale occurs at k ~ 3Kp. 

3.6. One dimensional structure functions of the velocity 



We have measured the structure functions of the ve- 
locity field to order p, 

where the velocity vector is decomposed in parallel and 
perpendi cular componen t s rela tive to the displacement 
vector t. iShe fc Leveaud (|1994H . have determined based 
on very general assumptions that these functions should 
scale oc within the inertial range of fully developed 
incompressible turbulence, where 

( p = 1 + 2-2(1)^ (15) 

These predictions have been extended analytically 
to i ncompressible MHD turb ulence (IGrauer et al.1 
19941) and verified num erically (|Politano et al.l 119951 : 
Muller fc Biskamdl2000l ). The compressible MHD case 
has been studied numeric ally by iBoldvrev et al.l (|2002b) 
and iPadoan et al.l (|2004[) in the context of supersonic 
molecular cloud turbulence, finding remarkably good 
agreement with analytical predictions. 

Figure E0 shows both s((£) and S^{t). We find that 
over the range 21A < £ < 82A, the structure func- 
tions are described by a power law with index (2 = 
1.025 ± 0.005. The high precision of this fit does not 
rule out systematic errors, but it does imply that some 
degree of scale-invariant behavior has been resolved in 
our simulations. The range over which the fit is valid 
occurs at smaller s cales than what was reported by 
iKritsuk et all (I2007D (32 A to 256 A) bu t very similar to 
those found bv IBoldvrev et all (|2002bf ) (10A to 90A). 
The slope (2 = 1.025 is slightly s teeper than what 
was observed by IKritsuk et al. (2007), and considerably 
steeper than the v alue of 2/3 which follows from the 
iKolmogorovl (|1941f ) theory. We also observe that the ra- 
tio S2 I s\ = 2.5 is lar g er by a factor of 2 than what was 
seen in IKritsuk et al.l (|2007l ). whose result is in closer 
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0.05 0.10 0.15 0.20 

Fig. 9. — Shown is the two-dimensional, second order structure 
function of the velocity field Sj (^-Li^ll)' The offset vector t is 
decomposed into components parallel (x-axis) and perpendicular 
(y-axis) to the local magnetic field, ^ (B ( r + + B(r)). 



agreement with the (|Kolmogorovlll94l"l ) theory. This ob- 
servation warrants farther attention, since if the discrep- 
ancy is not numerical in nature, but a feature of relativis- 
tic MHD turbulence then it may provide important hints 
in see king a relativistic extension to the iShe fc Leveaud 
(1994) model. 

Using about 10 9 sample pairs, we have been able to 
compute the slopes of higher order structure functions 
through p = 11. Figure [7] demonstrates the quality of 
these fits for each order. We have used the same win- 
dow to compute the slope of the higher order structure 
functions as for the S% case. Figure [8] shows the slopes 
£ p for the transverse and longitudinal structure functions 
for each p, normalized b y Ca- A lso shown are the predic- 
tions of IShe fc Leveaud (119941) for incompr e ssible non- 
relativistic hydro turbule nce. iGrauer et al.l (|1994f ) and 
IMuller fc Biskampl (|200fJ ) for incompressible MHD tur- 
bulence. We find a remarkable agreement between the 
IGrauer et al.l (|1994f) prediction and the longitudinal ve- 
locity fluctuations, even at the highest order. The data 
for transverse velocity fluctu ations lie midway between 
the IMuller fc Biskampl (|2000h prediction. 

3.7. Scale- dependent anisotropy of the velocity field 

MHD turbulence in the presence of a strong back- 
ground field, known as Alfvenic turbulence, consists 
of an energy cascade mediated by interacting MHD 
waves. This situation is treat able with the forma l- 
ism of wave turbulence (see e.g. iZakharov et al.l ll992). 
whereby the resonant interactions between the MHD 
"free particles" (solutions to the linearized equations) 
are treated perturbatively. For incompressible MHD, the 
only free particles are the shear and pseudo Alfven waves. 
Analysis of the resonant n onlinear interaction o f these 
modes lead to the model of IGoldreich fc Sridharl (1995) 
which predicts a Kolmogorov-like energy spectrum oc 
^,-5/3 Their findings also included the so-called scale- 
dependent anisotropy with respect to the local mean 

field, /c|| oc fc 2 ^ 3 . This phenomenon is understood geomet- 
rically as the exaggerated distortion of eddies, becoming 
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Fig. 10. — The semi-major (in) and semi-minor (£j_) axes of the 
elliptical contours obtained from Figure[9] (circles), and the best fit 
line (solid) lu oc ^ 84 . For comparison we provide the prediction 

of IGoldreich fc Sridharl 019950 (\\ oc £ 2 ^ 3 (dashed) and a slope of 
unity (dash- dotted). The vertical line marks the scale I = L/10 at 
which the magnetic and kinetic energy are in equipartition. 

more apparent at smaller sca les. The numer i cal ver i- 
fication was first obtained by ICho fc Vishniad (|2000aft . 
who presented two methods of numerical measurement 
corresponding to the correct geometrical interpretation. 
Other numerical studies have obs erved the same scaling 
in supersonic MHD t urbulence (|Cho fc Lazarianl 120031 : 
iBeresnvak et al.ll2005f) and fo rce-free relativistic Alfvenic 
MHD turbulence (|Chol 120051 ) , the latter having receive d 
prior analytical treatment bv lThompson fc BTa cs (1998). 
Here we provide the corresponding measurement for 
super- Alfvenic, compressible relativistic MHD turbu- 
lence. 

We u se the same method proposed by ICho fc Vis hniac 
(|2000a|) to measure the scale-dependence of eddy distor- 
tions, which relies on the eccentricities of level-surfaces 
for the second order structure functions 



: (|v(r- 
:<|B(r- 



v(r)| 2 ) 
-B(r)| 2 ) 



(16) 
(17) 



where t is decomposed into cylindrical coordinates ori- 
ented along the local mean field, which is defined for each 
pair of points as tj(B (r + £) + B(r)). 

In ex amining the shape of eddies, ICho fc Vishniad 
(2000a! used the structure functions for the velocity and 
magnetic field 5| and S 1 ^, finding slopes which were 
on average slightly larger than 2/3 for S% and slightly 
smaller for 5^. As depicted in Figure 1(b) coherent 
magnetic field structures exist out to only about 1/10 
of the domain size, meaning that provides a rather 
narrow window over which to measure the scaling. How- 
ever, the velocity field is coherent out the driving scale 
which affords a broad range of scales over which S*! scales 
robustly. 

Indeed, we observe excellent scaling in the structure 
functions of velocity, indicating that the slope £\\ oc i ^ 84 
is very robust. In fact, the validity of the fit is valid 
at least between 12 A and 300 A, the limiting factor in 
extending the scaling being due to data collection tech- 
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niques. The fact that the slope is steeper than the 2/3 
means that the eddy distor tions depend upon the scale 
more weakly than in the iGoldreich fc Sridharl (|1995H 
model. This may suggest corrections to the cascade dy- 
namics to account for rclativistic effects, but we forgo 
any conclusions on this matter until a detailed compar- 
ative study with the equivalent non-relativistic case has 
been completed. 

It is also very interesting that the same power law 
slope of 0.84 holds above and also below the equipar- 
tition scale. Under non - relativi s tic super- Alfvenic condi- 
tions, ICho fc Vishriiacl (|2000bD : ICho fc Lazarian] (|2003h 
both find similar contour shapes to those shown in Figure 
[9j However, in those studies the contour intercepts were 
not provided, meaning that the precise scaling behavior 
abo ve and below the equip artition scale is uncertain to 
us. iBeresnvak et al.l (|2005[ ) does provide these scalings 
for trans- Alfvenic, supersonic conditions, but only be- 
low the equipartition scale. We believe that a rigorous 
study of scaling above and below the equipartition scale 
in super- Alfvenic turbulence is required in order to iso- 
late the effects of compressibility, substantial magnetic 
field curvature, and also relativistic effects. 

4. CONCLUSIONS 

We have measured spectral and scaling properties of 
relativistically warm magnetohydrodynamic turbulence 
in the mildly compressible and super- Alfvenic regime. 
The numerical models were simulated at very high res- 
olution (1024 3 ) using Mara, a new second order Go- 
dunov code tuned for accurate and robust evolution of 
the RMHD equations in three dimensions. Our main pro- 
duction model was driven stochastically at large scales 
during quasi-steady evolution for about 6 light-crossing 
times of the domain. We find that: 

1. The magnetic energy is amplified from seed fields 
to 1.5% of the total fluid energy. The scale at which 
equipartition between magnetic and kinetic energy 
occurs is between 1/5 and 1/3 of the driving scale. 

2. At 1024 3 the power spectrum of velocity is dom- 
inated by a bottleneck, but not inconsistent with 
the Kolmogorov prediction of fc -5 / 3 . The power 
spectrum of density-weighted velocity p x ^ 3 v scales 
oc fc~ 5 / 3 over moderate wavenumb ers, consistent 
with the simple cascade model of iKritsuk et al.l 
ff2007h . 
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